#*************         NHANES       *************
#*************      Analyze Code    *************
#*************    线性回归模型讲解  *************
#*************                      *************
#### 0.准备好环境 ####
library(tidyverse)
library(gtsummary)
library(tidyr)
library(survey)

#### 1. 常规线性模型 ####
##### 1.1 生成演示用的数据信息 #####
# 生成 weight 变量，30 个男性 & 30 个女性
n = 30
weight.male <- rnorm(n, 75, 15)
weight.female <- rnorm(n, 55, 10)
weight <- c(weight.male, weight.female)

# 生成性别变量，30 个男性 & 30 个女性
sex.male <- factor(sample(c('male'), n, TRUE))
sex.female <- factor(sample(c('female'), n, TRUE))
sex <- c(sex.male, sex.female)


# 如果要随机生成 60 个男性和女性
# sex <- factor(sample(c('male', 'female'), n, TRUE))

# 生成身高变量
height <- 1*weight + rnorm(n * 2, 0, 5)

##### 1.2 生成演示用的数据信息 #####
# 构建简单的线性回归模型
data <- data.frame(height = height, weight = weight, sex = sex)
lm_1 <- glm(weight ~ sex + height, data = data) # Y ~ X1+X2, 英文波浪号
lm_2 <- glm(height ~ sex + weight, data = data) # Y ~ X1+X2, 英文波浪号
lm <- glm(height ~ sex + weight, data = data) # Y ~ X1+X2, 英文波浪号

summary(lm)
coefficients(lm)
confint(lm)
fitted(lm)
AIC(lm_1)
AIC(lm_2)
predict(lm)

##### 1.2 使用 gtsummary 优化输出 #####
lm_tbl <- tbl_regression(lm)

#### 2. 加权线性模型 ####
##### 2.1 提取数据模块 #####
demo.i <- read_xpt("2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)

##### 2.2 提取变量  & 去掉 NA 的行 & 衍生变量 #####
analyze.sample.data <- demo.i[,c('SEQN',
                                 'RIDRETH3', "RIDAGEYR", "RIAGENDR", "DMDEDUC2", "INDFMPIR", 
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]
# 一键去掉 NA 的行
analyze.sample.data.drop.na <- drop_na(analyze.sample.data) 
dim(analyze.sample.data.drop.na) 

# View(analyze.sample.data)

##### 2.3 生成复杂抽样 NHANES_design #####
NHANES_design <- svydesign(
  data = analyze.sample.data.drop.na, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~WTINT2YR
)

##### 2.4 线性回顾模型讲解 #####
# 年龄、教育和贫困指数之间的相关性
lm <- svyglm(INDFMPIR ~ RIDAGEYR + factor(DMDEDUC2), design = NHANES_design)

summary(lm)
coefficients(lm)
confint(lm)
fitted(lm)
AIC(lm)
predict(lm)

##### 2.5 使用 gtsummary 优化输出 #####
tbl_regression(lm)

# 显示更多的参数
tbl_regression(lm)%>% 
  add_global_p() %>%
  add_glance_table(
    include = c(nobs, AIC, BIC))

#### 3. 课后作业 ####
# 使用非加权的方式对 年龄、教育和贫困指数之间的相关性进行分析，并与加权结果进行对比

